Physiological Reports 



Physiological Reports ISSN 2051-817X 

ORIGINAL RESEARCH 

Combined reflectance spectroscopy and stochastic 
modeling approach for noninvasive hemoglobin 
determination via palpebral conjunctiva 

Oleg Kim 1 , John McMurdy 2 , Gregory Jay 3 , Collin Lines 4 , Gregory Crawford 4 & Mark Alber 1,4,5 

1 Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, Indiana, 46556 

2 Division of Engineering, Brown University, Providence, Rhode Island, 02912 

3 Department of Emergency Medicine and Division of Engineering, Brown University, Providence, Rhode Island, 02912 

4 Department of Physics, University of Notre Dame, Notre Dame, Indiana, 46556 

5 Department of Medicine, Indiana University School of Medicince, Indianapolis, Indiana, 46202 



Keywords 

Blood, hemoglobin, modeling, reflectance 
spectroscopy. 

Correspondence 

Mark Alber, Department of Applied and 
Computational Mathematics and Statistics, 
University of Notre Dame, Notre Dame, IN 
46556. 

Tel: 15746318371 
Fax: 15746314822 
E-mail: malber@nd.edu 

Funding Information 

This research was partially supported by the 
Gerber Foundation, National Science 
Foundation (NSF) grant DMS-0800612, and 
Rhode Island Science and Technology 
Advisory Council Grant. 



Abstract 

A combination of stochastic photon propagation model in a multilayered 
human eyelid tissue and reflectance spectroscopy was used to study palpebral 
conjunctiva spectral reflectance for hemoglobin (Hgb) determination. The 
developed model is the first biologically relevant model of eyelid tissue, which 
was shown to provide very good approximation to the measured spectra. 
Tissue optical parameters were defined using previous histological and 
microscopy studies of a human eyelid. After calibration of the model parame- 
ters the responses of reflectance spectra to Hgb level and blood oxygenation 
variations were calculated. The stimulated reflectance spectra in adults with 
normal and low Hgb levels agreed well with experimental data for Hgb con- 
centrations from 8.1 to 16.7 g/dL. The extracted Hgb levels were compared 
with in vitro Hgb measurements. The root mean square error of cross-valida- 
tion was 1.64 g/dL. The method was shown to provide 86% sensitivity esti- 
mates for clinically diagnosed anemia cases. A combination of the model with 
spectroscopy measurements provides a new tool for noninvasive study of 
human conjunctiva to aid in diagnosing blood disorders such as anemia. 
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Introduction 

Anemia is a condition in which the number of red blood 
cells or the amount of hemoglobin (Hgb) is below 
normal. According to World Health Organization, per- 
sons with concentrations less than 13 g/dL for men and 
12 g/dL for women are defined to be anemic. Anemia 
affects almost 3.5 million Americans, with millions more 
being undiagnosed, making it the most common blood 
disorder in the United States (Goodnough and Nissenson 
2004). Anemia can be induced by malaria, a mosquito- 



borne infectious disease caused by the Plasmodium para- 
site. Each year severe malarial anemia causes between 
190,000 and 974,000 deaths among children less than 
5 years (Murphy and Breman 2001) of age. Anemia is 
also the most common complication seen in 80% people 
with AIDS/HIV (Volberding et aL 2004). 

There exist different types of anemia ranging from mild 
and easily treatable iron and vitamin deficiency anemias 
to serious life-threatening aplastic and sickle cell anemias. 
The main anemia symptoms affecting physical function 
are fatigue, weakness, and fainting. Meanwhile, mild 
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anemia may occur without symptoms and can be detected 
only during a medical exam including a blood test. 

Current clinical methods to diagnose anemia rely on 
invasive determination of blood Hgb. Generally, the 
complete blood count (CBC) test is performed, which 
requires a venipuncture and a hematological analysis of 
the collected blood. Although the CBC test is very accu- 
rate, it has certain weaknesses such as a high cost of labor 
and instrumentation, and time-expensive analysis of the 
test, which frequently eliminates it from routine physical 
exams. In addition, because of the invasiveness of the test, 
the phlebotomist is exposed to risk of interacting with 
blood-born diseases. It is clear that a noninvasive method 
allowing fast and easy determination of Hgb is needed. 

There exist a number of noninvasive optical methods 
to indirectly measure the amount of Hgb and blood 
oxygenation levels in human tissue. These include retinal 
imaging (Rice et al. 2002), blood oxygenation monitoring 
using continuous NIR spectrometry (Liu et al. 1999; 
Benni 2002), photoplethysmography (Aldrich et al. 2002), 
ultraviolet and visible spectroscopies (Mourant et al. 
2005; Bigio and Mourant 1999), as well as fluorescence 
spectroscopy of oral tissue (Pavlova et al. 2009). 
(A detailed review of different noninvasive methods 
of total Hgb determination is presented in [McMurdy 
et al. 2008]). Tissue reflectance spectroscopy allows one 
to obtain optical properties of the tissue by measuring 
its reflectance spectra. Reflectance spectroscopy was 
previously implemented to characterize cutaneous tissues 
(Bigio et al. 2000; Mourant et al. 2005; Randeberg et al. 
2010, 2011; Kim et al. 2012; Klein et al. 2012), which are 
easily accessible, but present difficulties when monitoring 
Hgb signal due to interracial/interpatient melanin fluctua- 
tions, high scattering in turbid media, and weak signals 
due to skin tissue absorption at visible and near infra-red 
wavelengths. These factors can potentially influence the 
necessary instrumentation and inhibit the acceptance of 
noninvasive technology by clinicians. 

Human palpebral conjunctiva represents an easily 
accessible site for noninvasive Hgb level monitoring due 
to proximity of blood capillaries to an eyelid internal 
surface, which does not contain strong light absorbing 
chromophores such as melanin. 

There have been several studies conducted to evaluate 
Hgb level from conjunctiva spectra (Ernsting et al. 2001; 
McMurdy et al. 2006; Suner et al. 2007). The feasibility 
of applying image analysis algorithms to digital 
photographs of the conjunctiva to predict Hgb concentra- 
tions was examined with moderately successful results 
(Ernsting et al. 2001). Conjunctival reflectance spectra 
coupled with a partial least-squares (PLS) analysis and a 
discrete spectral region model were used (McMurdy 
et al. 2006) to derive Hgb concentration and compared 
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to in vitro measurements. Although these studies pro- 
vided some correlation between Hgb concentration and 
reflectance, they did not use a physiologically based 
model to reproduce the spectra of conjunctiva, and thus 
miss information about conjunctiva structural features 
which is important for spectral analysis when assessing 
Hgb levels. 

Because the tissue structure is complex, determining 
chromophore concentrations from the spectra involves 
the solution of an inverse problem, which is a nontrivial 
task (Pavlova et al. 2009; Phelps et al. 2010). The goal 
of this study was to describe a novel method combin- 
ing reflectance spectroscopy and a stochastic modeling 
approach of light propagation in multilayered tissue 
(Wang et al. 1995; Churmakov et al. 2002, 2003; 
Meglinski and Matcher 2002, 2003; Kim et al. 2012) to 
determine the reflectance of the palpebral conjunctiva. 
This will aid in evaluating the concentration of blood 
Hgb level and thus, in making a correct patient-specific 
medical diagnosis. 

Common approaches to model light propagation in 
biological tissues are based on the diffusion approxima- 
tion (DA) of a radiative transport equation and the 
stochastic Monte Carlo (MC) method. In comparison to 
MC method, DA is computationally less expensive and 
can be implemented for complex geometries. However, 
DA underestimates the light energy deposition in the epi- 
dermis and vessels when compared with MC (Zhang 
et al., 2005b). The difference is nonlinearly dependent on 
the wavelength, dermal blood volume fraction, vessel size, 
and depth. It was also shown (Flock et al. 1989) that 
the predictions made using DA became increasingly 
inaccurate as the albedo tended to zero and as the scatter- 
ing anisotropy factor became large (g>0.1) (Tuchin 
2007), which is typical for biological tissues where 
g~ 0.6—0.9. Meanwhile, MC method does not have these 
limitations. Another chief advantage of MC approach is 
that complex geometries and inhomogeneities can be 
modeled in a straight forward manner, and that a variety 
of physical quantities can be scored during the same run. 

The novelty of our extension of the MC-based 
modeling approach used in this study lies in its modified 
representation of the conjunctiva physiology (including 
spectral characteristics of conjunctival epithelium). Our 
novel solution method of the inverse problem quantifies 
Hgb level by minimizing the reflectance error over the 
broad wavelength range from 400 to 660 nm, whereas in 
previous models only several wavelengths were considered 
(Phelps et al. 2010). To study interaction of light photons 
with an eyelid, we developed an eye tissue model 
consisting of seven layers. In this model, each layer is 
characterized by a number of different parameters 
including chromophore concentrations, layer thickness, 
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refraction and anisotropy indices, and scattering and 
absorption coefficients. A stochastic Monte Carlo-based 
method is used to simulate light propagation through the 
tissue and conjunctival reflectance spectra are calculated 
and compared with experimental data. By fitting the 
measured and calculated reflectance curves, blood Hgb 
level is determined. 

The rest of the study is organized as follows. Materials 
and Methods section describes experimental procedure of 
human conjunctiva reflectance spectroscopy. Then, the 
Tarsal Conjunctiva Model provides details of the seven- 
layer tissue model. Next, the MC Method of photon 
propagation is described. The Results and Discussion 
section describes model simulation results, comparison 
between the simulated and experimentally measured 
spectra for different Hgb levels and discusses limitations 
and further extension of the approach. Finally, 
conclusions about the promising performance of the sug- 
gested noninvasive method and anemia sensitivity screen- 
ing are made. This justifies the practical application of 
the approach that can potentially result in a valuable 
diagnostic tool for medicine. 

Material and Methods 

To calibrate the model, we used the spectral data 
collected from prior patient data. Patients enrollment 
and collection of conjunctiva spectra are described in 
detail in our previous work (McMurdy et al. 2006). 
Briefly, 32 patients were enrolled at the Emergency 
Department at Rhode Island Hospital Providence, Rhode 
Island. Patients with a variety of complaints were 
enrolled to simulate a random patient population. 
Patients with acute cardiac symptoms, musculoskeletal 
injuries inhibiting comfortable relocation to the research 
location, symptoms of stroke, patients not yet seen or 



cleared by a physician, and critically ill patients were 
excluded from this study. Of the 32 consented patients, 
30 were able to participate (14 male and 16 female). 
Two patients were unable to participate; one due to 
nausea and abdominal pain during the enrollment 
procedure, and one because the conjunctiva could not 
be sufficiently exposed for spectroscopic measurements. 
Twenty-one of these patients were Caucasian, six were 
Hispanic, and three were African-American. Additionally, 
only patients with a minimum oxygen saturation of 90% 
were enrolled to exclude hypoxic patients during preli- 
minary evaluation of the efficacy of this technique. The 
mean oxygen saturation was 98.3 ± 1.5%. From CBC 
tests of 30 patients, the mean Hgb concentration was 
determined to be 13.9 ± 2.1 g/dL and the mean 
hematocrit was 40 ± 6.0. The correlation between Hgb 
and hematocrit values is shown in Figure 1A which is 
used as an input parameter for the model. Patient's con- 
junctiva was illuminated at a separation of ~15 cm from 
source element to conjunctiva. A quartz-tungsten-halo- 
gen source (Lot-Oriel, Surrey, UK) with total broadband 
power of 30 W was used for all measurements. The 
power level was chosen such that the procedure would 
not cause discomfort to a patient. The irradiance of this 
source did not exceed 0.1 W/(nm m 2 ) over the visible 
wavelength range. According to ANSI Z136.1 standards 
for incoherent radiation exposure, the procedure did not 
pose a radiation exposure risk (Charschan and Rockwell 
1999). All reflectance spectra (Fig. IB) were collected 
with an aperture size of 1-deg field of view (an aperture 
size of ~4.5 mm 2 ). The separation distance from 
spectrometer to conjunctiva was 7 cm. Although a 
smaller aperture would be more suited to collect reflec- 
tance spectra from single capillary vessels, the random 
movements of patients during acquisition did not permit 
this procedure. 




Hemoglobin, [g/dL] Wavelength, k [nm] 

Figure 1. (A) The relationship between measured hematocrit and hemoglobin values. Abscissa shows hemoglobin values in grams per decaliter 
and ordinate shows the corresponding hematocrit values in terms of the volume percentage. Symbols are data and the line is a linear 
regression. (B) Experimental reflectance spectra measured by a grating spectrometer in different patients with hemoglobin concentration 
between 7.3 and 16.7 g/dL. 
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Tarsal conjunctiva model 

Physiology of tarsal conjunctiva 

The conjunctiva is a mucosal membrane lining the eyelid 
(palpebral conjunctiva), the fold in the eyelid (forniceal), 
and the sclera (bulbar conjunctiva). Palpebral conjunctiva 
is an attractive location to diagnose anemia as it repre- 
sents a highly vascular area with a number of surface cap- 
illaries and the mucous membrane of the conjunctiva is 
highly transparent. Physiologically, palpebral conjunctiva 
is subdivided into margin conjunctiva, tarsal conjunctiva, 
and orbital conjunctiva. Tarsal conjunctiva overlying the 
internal surface of the tarsal plate is the best location for 
determining Hgb through spectroscopy measurements. 
This region has a thin epithelium layer containing no 
melanin, which does not interfere the extraction of Hgb 
signal from reflectance spectra. 

Figure 2 shows the histological image of a human low 
eyelid at tarsal plate location. It can be subdivided into 
several regions: conjunctival epithelium, tarsal plate, 
orbicularis oculi, subcutaneous tissue, dermis, epidermis, 
and stratum corneum on the outside of the eyelid tissue. 
Tarsal conjunctiva epithelium layer contains 2-3 layers of 
cylindrical cells and cuboidal cells. Tarsal plate contains 
collagen type I with rows of fibroblasts, which are fiber- 
forming cells in between collagen fibers. The orbicularis 
oculi is composed of striated muscle cells. It is responsible 
for blinking and squeezing eyelids shut. The eyelid skin is 
the thinnest in the body. It is represented by fat tissue, 
dermis, epidermis, and stratum corneum. Fat tissue con- 
sists of adipocytes, nerves, and blood vessels. Dermis 
mostly consists of supporting matrix which attract and 



retain water. Other components also include elastic and 
collagen fibers, blood vessels, hair follicles, sweat glands, 
macrophages, nerve endings, and lymphatic elements. The 
epidermis is comprised of several layers of keratinocytes 
with a superficial stratum corneum layer, which help to 
keep the skin hydrated by preventing evaporation of 
water. Conjunctiva derives its blood supply from the 
ascending branch of the posterior conjunctival artery. The 
major blood vessels of palpebral conjunctiva are located 
in the junction of tarsal conjunctiva and orbicularis oculi. 
There are also arterioles and capillary networks penetrat- 
ing orbicularis oculi, tarsal plate, subcutaneous tissue, and 
dermis. 

Parametrization of tarsal conjunctiva 

The structural complexity of a heterogeneous tissue is a fac- 
tor affecting light propagation (Jacques 1996; Verkruysse 
et al. 1999; Meglinski and Matcher 2002; Tuchin 2007). 
The spatial distribution of tarsal conjunctiva chromophore 
components including blood is not uniform (Ryan 1983; 
Renkin and Michel 1984). However, one can distinguish 
layers in which chromophore content is roughly constant. 
Therefore, as a first-order approximation, tarsal conjunc- 
tiva can be represented as a multilayered structure, with tis- 
sue components uniformly distributed within each layer. 
Some variability in thickness is expected from region to 
region of the conjunctiva and between individuals, and 
histological evidence in anecdotal human studies suggests 
this can be of order 30-70%. Nevertheless, for our study, 
we used these values for layer thickness to be typical for 
adults. 




A B 

Layer N |EM A DT 
1 =M= 
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Figure 2. (A) A multilayered tissue model approximates a lower eyelid with seven layers which optical properties are specified in Table 1 . 
Photons generated by a light emitter (EM) propagate in a tissue and the reflectance signal is measured by a spectroscope detector (DT). (B) 
Lower eyelid histology (http://histology.med.umich.edu/medical/eye) showing layers of conjunctival epithelium (1), tarsal plate (2), orbicularis 
oculi (3), subcutaneous fat (4), dermis (5), living epidermis (6), and stratum corneum (7). (C) The photograph of the data collection procedure 
using grating spectrometer. 
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We used a seven layers tissue model to simulate light 
propagation through palpebral tarsal conjunctiva with the 
model parameters given in Table 1. The absorption coeffi- 
cients for each layer takes into account the distribution of 
blood vessels, water and melanin and can be defined as 
(Meglinski and Matcher 2002) 

where \i l a and C; are the absorption coefficient and vol- 
ume fraction of the zth chromophore, N is the total num- 
ber of chromophores in the layer, and /i" is the 
absorption coefficient of Hgb- and water-free tissue. 

Experimentally measured absorption and scattering 
coefficients of the conjunctival epithelium (Nemati et al. 
1997, 1996) are approximated by least-squares polynomial 
fits as 

/£ E = AiX A + A 1 X i + A 3 ?. 2 + A 4 X + A 5 (2) 

where A x = 5.76328 x 10~ 13 A 2 = -5.05681 x 10~ 9 , A 3 = 
1.09816 x 10~ 5 A A = -0.00905007, A 5 = 2.607, and 

fif = Bxi 5 + B 2 X 4 + B 3 X 3 + B 4 X 2 + B 5 X + B 6 (3) 

where B; = 2.63396 x 10~ 13 , B 2 = -2.81535 x 10~ 10 , 
B 3 = -1.10282 x 10~ 6 , B 4 = 0.00236644, B 5 = -1.62947, 
and B 6 = 392.569. Here, the absorption and scattering fit 
errors, _R 2 , are 0.98 and 0.9, respectively. The units of fi^ E 
and /.i$ E are inverse millimeters and X is in nm. Polyno- 
mial coefficients A k and B k are calculated by minimizing 
the error, R 2 , of a polynomial fit, /i flt = J2k=o a kX k - To 
find a k , the equation for a fit in matrix notation, 
Mfi, = Aa, is solved for a given set of n experimental data 
points (/.,-, j.tj) and the solution vector a is calculated as, 
a = (A r A r 1 A r M fit (Nemati et al. 1997, 1996). 



Absorption spectra of Hgb, water and melanin are 
taken from the literature (Hale and Querry 1973; Wang 
et al. 1995; Prahl 1999). Absorption coefficient of Hgb 
/(J^ b (1/cm) is calculated as follows, 

/if (1) = 2.303 x e m (X) x x/M Hb (4) 

where e [L/(cm x mol)] is a Hgb molar extinction coeffi- 
cient (Prahl et al. 1989), x [g/L] is the number of Hb 
grams per liter, and M Hb = 64,500 (g/mol) is the Hgb 
gram molecular weight. The melanin absorption coeffi- 
cient is given by (Jacques and McAuliffe 1991), 

/(^ b e = 1.7 x lO'^r^Vm- 1 ] (5) 

and the absorption coefficient of the Hgb- and water-free 
tissue is approximated by (Jacques 1996), 

= 0.0244 + 8.53 exp[-(l - 154)/66.2] [mm -1 ] (6) 

with X given in nm. For blood containing layers, equa- 
tion (1) yields 

M«W = fC(X)7 Cbiocd(l ~S) + ^ 02 (X)yC bl00 dS 

+ /i H2 oW c H2 o(i - rc b i 00d ) + n° a (X) (7) 

(1 — yCbiood)(l — Ch 2 o) 

where y is the fraction of Hgb in blood, Cbiood is the 
blood volume fraction, and S is the oxygenation parame- 
ter. Assuming that Hgb is contained in erythrocytes only, 
y = <p H b Ht, where tp Hb is the volume fraction of Hgb in 
a single erythrocytes and Ht is the hematocrit level. To 
determine the hematocrit-Hgb relation, the linear regres- 
sion analysis based on 30 patients' CBC tests was used, 
which yielded the following Ht(Hgb) dependence, (see 
Fig. 1A) 

Ht[%] = 2.5829 x Hgb[g/dL] + 4.73[%] (8) 



Table 1. Optical and physiological parameters used for the tarsal conjunctiva model (Meglinski and Matcher 2002; Kakizaki et al. 2009; Cho 
et al. 2011; Kim et al. 2012). 



Layer 


Thickness, h jim 


9 


n 






Ch 2 o 


Cbiood 


S 




Ht, equation (8) 


CE 


20-70 


0.90 


1.34 


Equation (3) 














TP 


100-1500 


0.90 


1.4 


5 


0 


0.3 


0.3-0.6 


0.20- 


-0.99 


0.2-0.5 


OO 


2000-2500 


0.95 


1.4 


5 


0 


0.3 


0.3-0.6 


0.20- 


-0.99 


0.2-0.5 


ST 


1000-1500 


0.75 


1.44 


5 


0 


0.7 


0.06 


0.20- 


-0.99 


0.2-0.5 


D 


30-40 


0.80 


1.38 


35 


0 


0.6 


0.3 


0.20- 


-0.99 


0.2-0.5 


LE 


10-20 


0.80 


1.4 


45 


0.01-0.1 


0.2 


0 


0.20- 


-0.99 


0.2-0.5 


SC 


5-7 


0.86 


1.5 


100 


0 


0.05 


0 









Here, C me |, C H! o and Cbiood are the volume concentrations of melanin, water and blood. S is oxygen saturation parameter, and Ht is hemato- 
crit. CE, conjunctival epithelium; TP, tarsal plate; OO, orbicularis oculi; ST, subcutaneous tissue; D, dermis; LE, living epidermis; SC, stratum 
corneum. 
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Thus, equation (7) combined with equations (2-6), 

and (8) expresses the absorption of the blood containing 

layers as a function of wave length, A, Hgb level, Hgb, 

oxygenation, S, and known absorption coefficients [i™, 
Hbo 2 h 2 o J o 

MC Method 

The MC model employed the tracking of photon packets 
in a multilayer tissue mimicking human palpebral con- 
junctiva. The detailed description of the MC algorithm 
can be found in (Kim et al. 2012), (Wang et al. 1995). 
Briefly, the photon packet starts at a skin surface where 
radiation enters the medium and terminates when the 
packet is absorbed or leaves the tissue (Fig. 3). Inside the 
tissue the photon propagates in a random manner chang- 
ing its direction of motion each time it is scattered or 
internally reflected. The mean free path of the photon is 
calculated as S = — lnc;i/(/i f ), where £i is a uniformly dis- 
tributed random number from 0 to 1, (i t = /i s + ji a , with 
/( s and ji a being the scattering and absorption coefficients 
of the medium. At time zero each photon is assigned a 
unit weight. At each propagation time step, i, this weight 
is attenuated according to local absorption and scattering 
properties of the layer, such that a new weight is given by 
Wi = Vf;_i /( s //.( f . If Vfj becomes smaller than a prescribed 
threshold value w,h, then the photon is terminated 
and the next photon packet enters the tissue medium. 
If Wi is still larger than w,h, a new movement direction 
of the photon is calculated with azimuthal and 
longitudinal angles calculated as cp = 2n£ 2 > and 8 = 



Incident photon 
packet 


w = l 

t 


w 

n 

Reflectance 




f 7* 


h 


i / 


V.= (iy^r 7 





























Figure 3. Schematic diagram of photon propagation in a 
multilayered model tissue. 
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arccos[0.5(l+g 2 - (1 - g 2 )\l - g + 2gC 3 y 2 )/g] (Hen- 
yey and Greenstein 1941), where g is the anisotropy 
factor, and £ 3 is a random number uniformly distributed 
from 0 to 1. A photon packet directional vector is 
calculated according to (Prahl et al. 1989). 

By accumulating the total photon weight at a detector 
location, the reflectance at a given wavelength is 
calculated as a ratio of the total reflected photon weight 
to a total weight of incident photons. By repeating the 
process for different wavelengths, the reflectance spectra is 
numerically obtained. In this study, the reflectance spectra 
from 400 to 700 nm were simulated at intervals of 2 nm. 
To obtain the stable solution, 10 5 — 10 s iterations were 
used. Input model parameters determining the optical 
properties of each layer for a given wavelength included: 
layer thickness, absorption and scattering coefficients (n a 
and /( s ), oxygenation parameter (S), refractive index (n), 
and anisotropy scattering parameter (g). g(X) ranges from 
— 1 to 1 corresponding to perfect backscattering and 
forward scattering. It is defined as an averaged cosine 
value of the phase function, which gives the intensity 
distribution of scattered light with respect to the 
deflection angle 6. 

Calibration and predictive simulations 

The combined approach consists of two parts: model 
calibration (Fig. 4) and predictive simulations (Fig. 5). 
During calibration, first, the input parameters including 
the Hgb level (measured using CBC) and scattering coef- 
ficients are specified for each layer. Next, the iterative 
procedure starts. Each step of the iteration consists of (1) 
changing simulation blood concentration and oxygenation 
levels and thicknesses of tissue layers; (2) calculating 
absorption properties according to equations (1-8); and 
(3) running forward MC simulations for the specified 
layer scattering and absorption parameters. After that, 
simulated and experimental spectra are compared and the 
sum of squared residuals (SSR) is calculated for each set 
of parameters. Finally the minimum of SSR is found and 
the corresponding Cy ood , h b and S are determined, where 
hi is the thickness of the z'-th layer. 

Predictive modeling is based on the iterative solution 
of an inverse problem using forward MC simulations. By 
using the calibrated input parameters and by running 
stochastic simulations over a specific range of Hgb and 
oxygenation values, the corresponding Hgb and S values 
are derived by fitting a simulated spectrum to a measured 
one. To find the best-fit spectrum, the minimum of SSR 
is calculated for each set of parameter values. 

To estimate the performance of a predictive model 
combined with spectroscopic measurements, correla- 
tion between predicted and measured with CBC Hgb 
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Input: g, n, C mQ 

C mei (Table 1), 
Hgb from CBC test 



Experimental 
spectra of conjunctiva 
with known Hgb 



Iterate 

ft', c!, , , i = 1..7 
' blood ' 

and S over a given 
range (Table 1) 



I 



Calculate (J^for 
all tissue layers 



Run Monte Carlo 
simulation to 
calculate reflectance 



Calculate SSR 
between simulated 
and experimental 
reflectance spectra 



Find spectrum with 

min SSR:use h \ 
Cbiood,' = l-7 and 5 



Figure 4. A flowchart of the main steps of model calibration: (1) 
the input parameters including measured Hgb level and scattering 
coefficients are specified for each layer, (2) the iterative procedure 
begins with specifying blood concentration and oxygenation levels 
and layer thicknesses, (3) absorption properties are calculated 
according to equations 1-8, (4) Monte Carlo simulations are run 
for specified parameters, (4) simulated and experimental spectra are 
compared and the sum of squared residuals (SSR) is calculated for 
each set of parameters, (S) minimum of SSR is found and the 
corresponding CL, d , h h and S are determined. 



Input: n s , g, n, C mQ 

C mel (Table 1), 

hi C* i = 1 7 
w ' c blood'' 1 - / 

from calibration 



Iterate: Hgb 
from 6.0 to 19.0 g/dL 
and S from 0.4 to 0.99 



Calculate n a for 
all tissue layers 



I 



Run MC to calculate 
reflectance spectrum 
for given set of 
parameters 



Measured 
reflectance spectrum 
of conjunctiva 



I 



Calculate SSR 
between simulated 
and experimental 
reflectance spectra 



I 



Find spectrum with 
min SSR and 
determine 
Hgb and S 



Figure 5. A flowchart of the main steps of predictive modeling: (1) 
the input parameters including those obtained from calibration are 
specified for each layer, (2) the iterative procedure begins with 
specifying Hgb level and oxygenation levels, (3) absorption 
properties are calculated according to equations 1-8, (4) Monte 
Carlo simulations are run for specified parameters, (4) simulated 
and experimental spectra are compared and the sum of squared 
residuals (SSR) is calculated for each set of parameters, (5) the 
minimum of SSR is found and the corresponding Hgb and S are 
determined. 



levels is assessed. Cross-validation is performed using the 
leave-one-out method (Boser et al. 1992). Each round of 
cross-validation uses 29 spectra as a training set to predict 
a Hgb concentration corresponding to a single spectrum. 
To reduce variability, multiple rounds of cross-validation 
are performed for 30 spectra, and the validation results 
are averaged over the rounds. 

Results and Discussion 

Sensitivity analysis 

To assess the effect of variations in eyelid layers thicknesses 
and chromophore concentrations on a palpebral conjunc- 
tiva reflectance spectrum, the sensitivity analysis was per- 
formed by calculating the response of the resultant 
reflectance over 400-660 nm wavelength range to changes 



in a specific conjunctiva model parameter. First, variations 
in conjunctiva reflectance were evaluated for different 
thicknesses of a particular layer, whereas other parameters 
being fixed. Our simulations showed (Fig. 6) that for a 
given set of conjunctiva parameter values, the reflectance 
was most sensitive to variations in the thicknesses of the 
upper three layers (conjunctival epithelium, tarsal plate, 
and orbicularis oculi) resulted in up to 30% relative change 
in reflectance, whereas deeper layers contributed less than 
1% to the relative change in reflectance. Next, our 
simulations demonstrated that changes in blood volume 
fraction in tarsal plate and orbicularis oculi yielded 40-50% 
relative variation in conjunctiva reflectance (Fig. 7) and its 
spectrum is much less sensitive (<2%) to variations in 
blood volume fraction in subcuteneous tissue and dermis. 
The sensitivity of the conjunctiva spectrum to the blood 
oxygenation parameter S was evaluated by running simula- 



© 2014 The Authors. Physiological Reports published by Wiley Periodicals, Inc. on behalf of 
the American Physiological Society and The Physiological Society. 



2014 | Vol. 2 | Iss. 1 | e00192 
Page 7 



Reflectance Spectroscopy of Palpebral Conjunctiva 



0. Kim ef al. 




Figure 6. Palpebral conjunctiva reflectance spectra obtained with Monte Carlo simulations of photons transport in the tissue for various 
thicknesses of (A) conjunctival epithelium, h, = 20-70 ;im, (B) tarsal plate h 2 = 100-1000 ;im, (C) orbicularis oculi h 3 = 2000-2500 /im, (D) 
subcutaneous fat /i 4 = 250-1500 fim. Cy ood = 0.6,Cy ood = 0.3, S = 0.9, Hgb = 15.4 g/dL. The relative variation is shown in the inset figure. 




qI , , , , 1 — 
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Figure 7. Variations in Monte Carlo simulated palpebral conjunctiva 
plate layer Q,i OOC i = (0.08 - 0.3); (B) orbicularis oculi Cbi 00 d = (0.15 - 
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lectance spectra due to changes in blood volume fraction in (A) tarsal 
6). The relative variation is shown in the inset figure. 



tions for different S values varied from 0.3 to 0.9 (Fig. 8). 
As oxygenation increases in all tissue layers, Hgb becomes 
more oxygenated and the reflectance spectrum changes 
according to the total Hgb absorption. The corresponding 
changes are observed in the shape of the spectral curve: as 
the amount of oxygenated blood increases, the double peak 
at around 550 nm becomes more pronounced and the local 
peak at 440 nm shifts toward shorter wavelengths. Our 
simulations also showed that the total reflectance was not 
sensitive to changes of melanin in the epidermis of an eye- 
lid (Fig. 9). For typical skin layer thicknesses and chromo- 
phores concentrations (Table 1), relative reflectance 



variations were less than 1%, when C me i changed between 
0.01 and 0.1. Thus, the performed sensitivity analysis indi- 
cated that the structural complexity of the conjunctiva 
model can be partially reduced by neglecting factors which 
impact on the reflectance curve is relatively small. 

Model calibration and predictive 
simulations 

The extended MC model was calibrated using measured 
experimental spectra of a human conjunctiva. To achieve 
this, the model was iteratively run and model parameters 
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Figure 8. Monte Carlo simulations of palpebral conjunctiva 
reflectance spectra for different values of the oxygenation 
parameter, S. 

were varied to find the best fit to experimental data. For 
a known experimental spectrum, using the measured Hgb 
and hematocrit values, the model simulated spectra 
were calculated for various tarsal plate layer thicknesses 
(1—1.5 mm), blood oxygenation parameter S (0.2-0.9), 
and tarsal layer blood volume fraction (0.1-0.9). The cri- 
terion for the best fitted simulated spectrum was to mini- 
mize the SSR ^a[R s (a) - R m {).)] 2 , where R S (X) and 
R m (X) are the simulated and measured reflectance values 
at the wavelength of ).. Thus, the best fit parameters were 
determined and the mismatch was minimized. After cali- 
bration, the MC model was used to simulate reflectance 
spectra for various Hgb levels in the range from 400 to 
660 nm, which corresponds to Hgb level variations in a 
person whose conjunctiva spectrum was used for calibra- 
tion. 

Light attenuation in deep biological tissues depends on 
the effective attenuation coefficient defined as /i eff = 



^mel varia li£fts>' 




400 450 500 550 600 
Wavelength, X [nm] 



650 



Figure 9. Monte Carlo simulations of palpebral conjunctiva 
reflectance spectra for different values of melanin concentration 
C me | in the eyelid epidermis. Reflectance spectrum is not sensitive to 

Oriel- 



\/3^ a (X + /<), where n' s = ^(1 - g). To evaluate light 
attenuation of different conjunctiva layers, ^i eff was calcu- 
lated for different wavelengths using the calibrated input 
parameters and equations (1-8). Attenuation of different 
layers is shown in Figure 10 for Hgb values 8.1 and 
16.4 g/dL. As the Hgb concentration increases, the relative 
attenuation of blood containing layers also increases. 
Thus, for tarsal plate and orbicularis oculi, ^ e ff increases 
1.5 times at 1 = 800 nm and 3.8 times at X = 415 nm, as 
Hgb increase from 8.1 to 16.4 g/dL. 

Figure 1 1 A-D demonstrates the capability of the model 
to make predictions of the conjunctiva reflectance for a 
given set of tissue parameters and a particular Hgb level 
determined in single-spectrum calibration. Once the 
conjunctiva parameters were set (Hgb = 8.1 g/dL, h 1 = 




ill 10 



500 600 700 

Wavelength, X [nm] 



500 600 700 

Wavelength, X [nm] 



Figure 10. The effective coefficient of attenuation for different eyelid tissue layers calculated according to equation (1): conjunctival epithelium 
(1), tarsal plate (2), orbicularis oculi (3), subcutaneous fat (4), dermis (5), living epidermis (6), and stratum corneum (7). (A) S = 0.6, 
Hgb = 8.1 g/dL, cp Hb = 0.4; (B) 5 = 0.9, Hgb = 16.4, (p Hb = 0.4. 
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Figure 11. Monte Carlo simulation results for reflectance spectra obtained with the seven-layered model of a palpebral conjunctiva and 
comparison with experimental data for different levels of hemoglobin: 8.1 g/dL (A), 14.5 g/dL (B), 15.4 g/dL (C), and 16.7 g/dL (D). The 
results of model simulations are shown by circles and experimental spectra are shown by lines. Error bars show the standard deviation 
of experimentally measured spectra of an individual patient. Different colors correspond to different patients having the same hemoglobin 
level. 



20 /mi, h 2 = 100 /mi, h 3 = 2 mm, h 4 = 1 mm, h s = 300 
/mi, h 6 = 200 /mi), h 7 = 7 /mi, Cj; lood = 0.6, and Cj^ lood = 
0.3), the conjunctiva reflectance spectra were calculated 
for various Hgb levels up to Hgb = 16.7 g/dL and com- 
pared with experimental data (Fig. IB). Here, the mean 
experimental values are shown by lines and standard devi- 
ations are shown by error bars, based on three individual 
conjunctiva spectra, measured 30 sec apart. Variations in 
errorbars (Fig. 1A vs. Fig. ID) of experimental recordings 
can be attributed to the fact that averaging was performed 
over a small number of spectra collected from conjuncti- 
vas of both eyes of a patient. One can expect that 
patient's left and right eyes' conjunctivas vary in tissue 
structure, Hgb concentrations as well as blood content. 
The orientation of the light emitter and detector with 
respect to a conjunctiva surface can also slightly affect the 
measurements. Meanwhile, the model correctly predicts 
both qualitative and quantitative behavior of the mea- 
sured reflectance. The correlation coefficients between the 
modeled and measured spectra were found to be in the 
range from 0.992 to 0.996. The simulated and experimen- 
tal curves agreed within the experimental uncertainty. 
The SSR values were found to be 0.19, 0.3, 0.18, and 0.17 



for (Hgb = 8.1, 14.5, 15.4, and 16.7 g/dL) where each 
SSR value was calculated with respect to a spectrum aver- 
aged over the inter- and intrapatient spectra over the 
wavelengths from 400 to 660 nm. 

Figure 11B, C also provide the degree of the intra- and 
interpatient spectral variations for the same Hgb values, 
which were found to be over 10-20% for the wavelengths 
below 580 nm. Because individual spectra were taken 
30 sec apart, the higher intrapatient spectral variations in 
this experiment are most likely due to small shifts in 
sampling locations from one measurement to the next 
rather than a temporal change in the conjunctiva charac- 
teristics within this timeframe. Interpatient variations also 
points out the importance of collecting multiple spectra 
from multiple palpebral conjunctiva locations of different 
individuals to perform a more careful uncertainty analy- 
sis. These will be accounted for in future studies. 

To find the degree of correlation between the predicted 
and CBC measured Hgb concentrations, we followed the 
procedure mentioned in the previous subsection and 
shown in Figure 5. Model calibration was performed by 
using leave-one-out cross-validation with 29 spectra used 
as a training set to predict a Hgb concentration 
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corresponding to a single spectrum. For each spectrum of 
a specific patient, tissue parameters were determined and 
Hgb and oxygenation levels were varied to simulate 
multiple spectral curves. Hgb was then extracted by 
minimizing SSR value for each measured spectrum. This 
was repeated 30 times such that each spectrum was used 
once as the validation dataset. After that, the extracted 
Hgb values and Hgb measured in vitro were cross 
correlated. The Pearson correlation coefficient was found 
to be 0.61. The root mean squared error of cross-valida- 
tion (RMSECV) was 1.67 g/dL (Fig. 12). For the seven 
patients whose CBC test showed they were clinically 
anemic, our method predicted six patients to be clinically 
anemic yielding the sensitivity of 86%. The results of the 
model-based algorithm are also presented as Bland- 
Altman plot in Figure 13. 

From the model simulations, the oxygenation 
parameter was found to be equal to 0.89, implying that 
tarsal plate averaged blood oxygenation level was 89%. 
Although it was not possible to measure the blood 
saturation level directly in the eyelid tissue, it is 
worthwhile estimating it. Pulse oximetry measurements 
provided oxygen saturation of arterial blood (Sa0 2 ) to 
vary from 94 to 100% with the mean value of 98.3%. The 
saturation of venous blood (Sv0 2 ) varies from 72% to 
75% (Sheinberg et al. 1992). Assuming that in the tissue 
layer oxygenated and deoxygenated blood present in equal 
parts, the weighted average of the total oxygen saturation 
varies from 85% to 88% which estimates the value of S 
within 2-7% accuracy. 




10 12 14 16 

CBC Hemoglobin, [g/dL] 



Figure 12. Prediction plot of Monte Carlo (MC)-based model 
versus invasive complete blood count (CBC) derived hemoglobin 
concentration. The correlation coefficient is 0.61 and 
RMSECV = 1 .64 g/dL. RMSECV, root mean squared error of 
cross-validation. 
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Figure 13. Bland-Altman plot of the model predicted hemoglobin 
versus the complete blood count (CBC) derived hemoglobin 
concentration. Solid line is at the mean difference of CBC/Model 
predicted hemoglobin (-0.1 g/dL) and dashed lines are drawn at 
±2 SD (±4.39 g/dL) to represent limits of agreement. 



Discussion 

The goal of this study was to introduce a new method of 
noninvasive assessment of blood Hgb by combining 
reflectance spectroscopy of human conjunctiva and 
stochastic model of light propagation through the 
conjunctiva tissue. The method has been used to evaluate 
Hgb levels and anemia conditions and the calculated 
anemia sensitivity was found to be 86%. We should 
notice that the obtained anemia sensitivity measure was 
calculated using one negative result, which was due to the 
relatively small number of anemic patients (seven) and 
the overall recruited patient population (30). Larger 
population study will be needed to improve anemia 
sensitivity measure. 

Due to the presence of multiple parameters characteriz- 
ing the tissue and high computational cost of the MC 
simulations, the difficulty in using present MC simulation 
was in adjusting the simulated reflectance spectra to the 
measured ones. The total computational cost of finding 
solution of the inverse problem for one spectrum was 
~ 24 h using 3.2 GHz processor computer. This corre- 
sponds to 30 parameter variation simulations over 200 
wavelengths using 10 5 -10 6 photons for a single-reflectance 
spectrum. Although we did not notice that results were 
strongly affected by the type of the iteration procedure, a 
more detailed analysis is needed to optimize the choice of 
the iteration algorithm in the future. Meanwhile, although 
it is essential to minimize the SSR over the entire spec- 
trum to accurately fit the reflectance curves, focusing on 
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specific characteristic wavelengths such as local extremum 
or isosbestic points can facilitate predictive simulations by 
decreasing the number of simulated wavelengths. How- 
ever, a detailed comparative analysis of how the reduced 
spectrum affects Hgb calculations would be the goal of 
future studies. 

Variations in conjunctiva parameters from person to 
person and spatial and wavelength dependencies of 
optical properties are not well described in the literature. 
From a few anecdotal human studies, 30-70% conjunc- 
tiva thickness variations can be estimated, however, more 
data should be collected to provide a reliable statistics. 
Meanwhile, from our simulations, the variation in layers' 
thicknesses reached 30-100%. Additional studies of 
human eyelid tissue properties are required to specify the 
proper distribution of absorption and scattering values. 
Conjunctiva spectral measurements accompanied by CBC 
tests of a specific individual over a longer timeframe 
whose Hgb varies in time will permit to narrow down 
model parameter uncertainty. These can be achieved by 
performing measurements during surgery or postsurgery 
recovery when blood properties can change over time 
(Luz and Rodrigues 2004). Further parameter optimiza- 
tion can be probably achieved using neural network algo- 
rithms although there may be unexpected errors due to 
local minima (Zhang et al., 2005a). 

The computational cost can be reduced by using 
GPU-based parallel computing which can provide 
acceleration from 340 to 1000 times (Alerstam et al. 
2008; Doronin and Meglinski 2011). Another potential 
reduction in the computational costs could be achieved 
by using adequate simplification of the model based on 
the sensitivity analysis performed in this study. As we 
have shown, conjunctiva spectra are mostly affected by 
the parameters of conjunctiva epithelium, tarsal plate, 
and orbicularis oculi. This means that potentially, the 
optical model of the low eyelid tissue can be simplified 
by considering only the aforementioned three layers, 
which would be enough to simulate the reflectance spec- 
trum of the palpebral conjunctiva without loss of accu- 
racy. As we have shown, the multilayered planar tissue 
structure provides very good spectral approximation 
results. However, further extension of the model would 
be the use of a higher order spatial approximation of 
the layer boundaries, such as the one used in (Meglinski 
and Matcher 2002). Notice that the use of a simpler 
model such as Kubelka-Munk theory-based model (KM) 
(Kubelka 1948) describing radiation transfer in a diffuse 
scattering medium, is inappropriate in our case. The KM 
is formulated as a two-flux model which considers only 
forward and backward fluxes propagating in an infinite 
slab of a certain thickness. High-order approximations of 
KM has been also reported (Meador and Weaver 1979; 



Jacques and Prahl 1987; Yang et al. 2004), however, their 
usefulness is limited because all these models are gener- 
ally restricted to a simple slab geometries, diffuse irradi- 
ance and isotropic scattering, which is atypical for 
biological tissue optics problems. The relative simplicity 
and speed of KM models are achieved at the expense of 
accuracy, which requires a more detailed analysis of the 
structure and optical properties of different tissue layers. 
Thus, in our case, the use of a more detailed and 
complex model based on the inverse MC method and 
not having the aforementioned limitations, is more 
appropriate. 

Conclusions 

Hgb and blood oxygenation levels are frequently used as 
indicators of human conditions and different blood dis- 
eases such as anemia and hypoxemia. Complementing 
visual inspection of a patient's conjunctiva with a tool for 
spectroscopical analysis of the palpebral conjunctiva can 
potentially provide a reliable noninvasive method to diag- 
nose these conditions. In this study, we have extended 
our previous studies of conjunctiva reflectance spectros- 
copy by coupling it with a new modeling approach to 
predict spectra of patients' conjunctiva. 

We have described a method which uses a combination 
of reflectance spectrometry and stochastic simulation of 
photon movement in eyelid tissue to allow assessment of 
a human blood Hgb level. A new seven-layer tissue model 
based on histological and optical coherence tomography 
studies was developed to account for biological details of 
the human eyelid structure. Different variations in optical 
parameters were considered to assess the sensitivity of the 
model. In particular, we showed that the conjunctiva 
reflectance is sensitive to variations in the thickness of the 
conjunctiva epithelium, tarsal plate, and orbicularis oculi 
as well as to variations in blood volume fractions of the 
tarsal plate and orbicularis oculi. 

The experimental conjunctival spectra of adults includ- 
ing anemic patients were used to calibrate the model and 
run predictive simulations. We showed that the simulated 
spectra represented a good approximation to the experi- 
mentally measured spectra of a human palpebral conjunc- 
tiva over the range of different Hgb levels from 8.1 to 
16.7 g/dL. The combined predictive simulations and mea- 
sured spectra of 30 patients permitted us to determine 
Hgb levels and compare them with Hgb measured in 
vitro using CBC test. The results of the study demonstrate 
that combination of measurements using a portable spec- 
troscopic device and a fast computer-based model analysis 
software can potentially result in development of an 
effective novel diagnostic tool for planning the treatment 
of patients. 
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